load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"
load "/public/users/wrf/Codes/WRF_scripts/utils/substring.ncl"
load "/public/users/wrf/Codes/WRF_scripts/utils/WRF_dyy_contributed.ncl"

; open csv file to generate the polygon
csvfile = "/public/users/wrf/Data/Northeast_rectangle.csv"
polygon = readAsciiTable(csvfile,2,"float",1)
dim_polygon = dimsizes(polygon)

; open wrfout file to get the domain lat/lon information
wrfout_file = "/public/users/wrf/wrfout/2012111420/wrfout_d01_2012-11-14_12:00:00"
a = addfile(wrfout_file + ".nc", "r")
xlong = wrf_user_getvar(a, "XLONG", 0)
xlat  = wrf_user_getvar(a, "XLAT", 0)
dim_domain = dimsizes(xlong)
print(dim_domain)
lonlat_recs_pool = new((/252, 2/), "float")
llcorner_recs_pool = new((/252, 5, 2/), "float")
ij_recs_pool = new((/252, 2/), "integer")
rec_count = 0

; determine the max and min i,j
lon_bdry = (/min(polygon(:,0)), max(polygon(:,0))/)
lat_bdry = (/min(polygon(:,1)), max(polygon(:,1))/)
print(lon_bdry)
print(lat_bdry)

ijCorner1 = wrf_user_ll_to_ij(a, min(polygon(:,0)), min(polygon(:,1)), True)
ijCorner2 = wrf_user_ll_to_ij(a, max(polygon(:,0)), max(polygon(:,1)), True)

print(ijCorner1)
print(ijCorner2)

do i = ijCorner1(0), ijCorner2(0), 1
    do j = ijCorner1(1), ijCorner2(1), 1
        print("i= " + i + " j= " + j)
        if (is_point_in_polygon((/xlong(j,i), xlat(j,i)/), polygon, dim_polygon(0))) then
            print("point ("+tostring(xlong(j,i))+", "+ tostring(xlat(j,i)) + ") is in the Northeast Rectangle.")
            lonlat_recs_pool(rec_count, :) = (/xlong(j,i), xlat(j,i)/)
            ij_recs_pool(rec_count, :) = (/i, j/)
            llcorner_recs_pool(rec_count, 0, :) = wrf_user_ij_to_ll(a, i+1, j+1, True)
            llcorner_recs_pool(rec_count, 1, :) = wrf_user_ij_to_ll(a, i+1+0.5, j+1+0.5, True)
            llcorner_recs_pool(rec_count, 2, :) = wrf_user_ij_to_ll(a, i+1-0.5, j+1+0.5, True)
            llcorner_recs_pool(rec_count, 3, :) = wrf_user_ij_to_ll(a, i+1-0.5, j+1-0.5, True)
            llcorner_recs_pool(rec_count, 4, :) = wrf_user_ij_to_ll(a, i+1+0.5, j+1-0.5, True)
            rec_count = rec_count + 1
        end if
    end do
end do

lonlat_recs = lonlat_recs_pool(0:rec_count-1, :)
llcorner_recs = llcorner_recs_pool(0:rec_count-1, :, :)
ij_recs = ij_recs_pool(0:rec_count-1, :)

printArray = sprintf("%7.3f", lonlat_recs(:,0)) + ", " + sprintf("%7.3f", lonlat_recs(:,1)) \
    + ", " + sprinti("%4i", ij_recs(:,0)) + ", " + sprinti("%4i", ij_recs(:,1)) \
    + ", " + sprintf("%7.3f", llcorner_recs(:, 0, 0))  + ", " + sprintf("%7.3f", llcorner_recs(:, 0, 1)) \
    + ", " + sprintf("%7.3f", llcorner_recs(:, 1, 0))  + ", " + sprintf("%7.3f", llcorner_recs(:, 1, 1)) \
    + ", " + sprintf("%7.3f", llcorner_recs(:, 2, 0))  + ", " + sprintf("%7.3f", llcorner_recs(:, 2, 1)) \
    + ", " + sprintf("%7.3f", llcorner_recs(:, 3, 0))  + ", " + sprintf("%7.3f", llcorner_recs(:, 3, 1)) \
    + ", " + sprintf("%7.3f", llcorner_recs(:, 4, 0))  + ", " + sprintf("%7.3f", llcorner_recs(:, 4, 1))

output_filename = "/public/users/wrf/Data/Northeast_domaian_mask_d01.csv"
print(printArray)
asciiwrite(output_filename, printArray)

print(rec_count)
